clear;
clc;
tic;
%参数
N = 50;    
n = 100;
xk = linspace(0,1,N)';     %初始中心点
xe = linspace(0,1,n)';     %评估点

%%
F = @(x) sin(5*pi*x) .* exp(-2*x);
y=F(xe);

%% LOOCV优化epsilon
epsilon_list = linspace(0.01, 20, 100); % 候选epsilon范围
cv_errors = zeros(size(epsilon_list));
% 预先计算距离矩阵（n x N）
DM = pdist2(xe, xk); % 或者使用abs(xe - xk')，但需要确认维度
for i = 1:length(epsilon_list)
    epsilon = epsilon_list(i);
    Phi = exp(-epsilon * DM.^2); % 生成基函数矩阵
    
    % 计算系数
    w = pinv(Phi) * y;
    
    % 计算帽子矩阵的对角元素
    H = Phi * pinv(Phi);
    H_ii = diag(H);
    
    % LOOCV残差
    residual_loo = (y - Phi * w) ./ (1 - H_ii);
    cv_errors(i) = mean(residual_loo.^2);
end
% 选择最优epsilon
[~, idx] = min(cv_errors);
epsilon_opt = epsilon_list(idx);
fprintf('LOOCV优化结果：最优epsilon=%.4f\n', epsilon_opt);
% 使用最优epsilon重新生成Phi矩阵
epsilon = epsilon_opt;
Phi = exp(-epsilon * DM.^2);

rbf=@(r)exp(-epsilon*(r).^2);

%% 
%构建评估矩阵
dist_matrix=pdist2(xe,xk);
A=rbf(dist_matrix);

%使用k折交叉验证，选择lambda
[W,FitInfo]=lasso(A,y,'CV',5);   %使用5折交叉验证

%选择最优lambda
lambda_option=FitInfo.Lambda1SE;
fprintf('最优lambda为：%.4f\n',lambda_option);

%使用最优lambda重新训练
w=W(:,FitInfo.Index1SE);    %对应权重
select_idx=find(w~=0);      %选出不为零的权重
select_xk=xk(select_idx);   %通过权重的索引，找出要选择的中心点，即选择的基函数
  
%构造稀疏模型
A_sparse=A(:,select_idx);
w_sparse=A_sparse\y;        %最小二乘优化权重
fe=A_sparse*w_sparse;

%误差分析
rmse=sqrt(mean((y-fe).^2));         %均方根误差
error=max(abs(fe-y));        %最大绝对误差
fprintf('均方误差rmse:%.2e\n最大绝对误差:%.2e\n',rmse,error);
fprintf('使用基函数数量%d\n',length(select_xk));

%% 可视化

%  LOOCV误差曲线
figure;
semilogy(epsilon_list, cv_errors, 'b-', 'LineWidth', 1.5);
xlabel('形状参数γ'); ylabel('LOOCV误差');
title('LOOCV误差随γ变化曲线');
grid off;


figure;
plot(xe,y,'k-','LineWidth',1.5);
hold on;
plot(xe,fe,'r--','LineWidth',1.5);
xlabel('x');
ylabel('f(x)');
title('使用交叉验证的lasso稀疏化RBF');
legend('真实解','lasso稀疏解');
grid on;
toc;